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Adaptive radiations present fascinating opportunities for studying the evol- 
utionary process. Most cases come from isolated lakes or islands, where 
unoccupied ecological space is filled through novel adaptations. Here, we 
describe an unusual example of an adaptive radiation: symbiotic mussels 
that colonized island-like chemosynthetic environments such as hydrothermal 
vents, cold seeps and sunken organic substrates on the vast deep-sea floor. Our 
time-calibrated molecular phylogeny suggests that the group originated and 
acquired sulfur-oxidizing symbionts in the Late Cretaceous, possibly while 
inhabiting organic substrates and long before its major radiation in the 
Middle Eocene to Early Oligocene. The first appearance of intracellular and 
methanotrophic symbionts was detected only after this major radiation. 
Thus, contrary to expectations, the major radiation may have not been trig- 
gered by the evolution of novel types of symbioses. We hypothesize that 
environmental factors, such as increased habitat availability and/ or increased 
dispersal capabilities, sparked the radiation. Intracellular and methanotrophic 
symbionts were acquired in several independent lineages and marked the 
onset of a second wave of diversification at vents and seeps. Changes in habitat 
type resulted in adaptive trends in shell lengths (related to the availability 
of space and energy, and physiological trade-offs) and in the successive 
colonization of greater water depths. 
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1. Introduction 

Adaptive radiation is broadly defined as the rapid diversification of species and of 
their adaptations to the environment in response to natural selection and ecologi- 
cal opportunities [1]. The radiations of Darwin's finches and African Rift Lake 
cichlids are well known and have become popular beyond academic research, 
but these examples may just be the 'tip of an evolutionary iceberg 7 , because 
adaptive radiations can take place on a broad range of time scales and taxonomic 
levels, and may even include the Cambrian explosion of life [2]. The fauna living 
at hydrothermal vents and cold seeps in the deep sea represents another remark- 
able, yet poorly understood case of radiation. These habitats are characterized 
by extreme physico-chemical features and a scarcity of primary photosynthetic 
production, and the animals living there thrive because of their symbiotic 
relationships with bacteria able to use sulfide and/ or methane as energy sources 
[3]. Initially considered as archaic faunas that survived major extinction events [4], 
molecular age estimates and fossil records suggested that most of the modern 
vent and seep animals appeared during a short time interval between the Late 
Mesozoic and the Early Cenozoic [5,6]. 
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Mussels of the bivalve family Mytilidae are particularly suit- 
able as model organisms to study the roles of ecological 
opportunities, symbioses and other adaptations in the evolution 
of deep-sea chemosymbiotic faunas. Indeed, they dominate 
many vents and seeps, but are also common in other sulfide- 
rich habitats such as whale carcasses and sunken wood [7], 
which are considered to be evolutionary stepping stones to 
deep-sea vents [8]. In addition to this diversity of habitats, 
deep-sea mussels have a remarkable range of symbiotic types, 
including intracellular and extracellular symbionts, and the abil- 
ity to host either or both sulfur-oxidizing and methanotrophic 
symbionts [7]. 

We investigated the evolutionary history of deep-sea sym- 
biotic mussels by estimating phylogenetic relationships from 
a comprehensive dataset including all known lineages and 
80% of the known species, collected from virtually all ocean 
basins. In an analysis of five gene fragments, including mito- 
chondrial and nuclear DNA, calibrated with three reliable 
fossils, we (i) characterized speciation rates through time; 
(ii) reconstructed the evolution of habitat use (environmental 
type and depth), body size, symbiont type (sulfur and /or 
methane oxidizer) and the degree of the physiological inte- 
gration of the symbiont with the host (intracellular versus 
extracellular symbioses); and (iii) evaluated the impact of 
these biological and ecological factors on speciation rates, 
and their respective roles in the evolution of deep-sea 
symbiotic mussels. 

2. Material and methods 

(a) Sampling 

Fourteen species of vent and seep mussels from the western 
Pacific were sampled using the manned submersible Shinkai 
6500 and remotely operated Hyper Dolphin vehicle (see elec- 
tronic supplementary material, table SI). From experimental 
bone deployments in Japanese waters, we collected the three 
undescribed species Idas sp SAL4, Idas evolutionary significant 
unit (ESU) D and Idas ESU R [9]. Upon recovery, pieces of the 
gills of Bathymodiolius aduloides and the three undescribed species 
were fixed with 2.5% glutaraldehyde in filtered seawater for 24 h 
and preserved in filtered seawater with 10 mM sodium azide at 
4°C. Remaining tissues and other specimens were fixed in 99% 
ethanol for DNA analysis. We also analysed alcohol-fixed foot 
tissues of 15 species associated with sunken organic substrates 
collected in the western Pacific during cruises of the Tropical 
Deep-sea Benthos program [9,10]. Finally, we included alcohol- 
fixed specimens from the Atlantic Ocean and Mediterranean 
Sea, including two species from organic falls, six species from 
hydrothermal vents, five species from cold seeps and the subtidal 
outgroup species Modiolus modiolus (Linnaeus, 1758). 

(b) Molecular analyses 

Template DNA from feet and gills were extracted to analyse host 
and symbiont genes, respectively, using the QIAamp DNA Micro 
Kit (Qiagen, USA). A fragment of the small subunit 16S rRNA 
was amplified to characterize symbionts in 10 species that had not 
been studied previously (see electronic supplementary material, 
table SI). Fragments of mitochondrial COI, NADH4 and 16S, and 
nuclear 28S and histone 3, were amplified from host species for 
phylogenetic analysis. Polymerase chain reactions (PCRs) were per- 
formed using the Ex Taq PCR kit (TaKaRa, Japan). Forward and 
reverse primers (0.2 |ulM each; see electronic supplementary 
material, table S2) and less than 1 juug of DNA template were 
added to reaction mixtures. PCR products were generated by an 



initial denaturing step of 4 min at 94°C followed by 35 cycles at 
94°C for 1 min, 55°C (50°C for COI and NADH 4 ) for 1 min and 
1 min at 72°C, and by a final extension at 72°C for 7 min. They 
were purified using Wizard SV Gels and the PCR Clean-Up 
System (Promega, Madison, WI). PCR products of symbiont 16S 
were cloned into the pCR-TOPO vectors using a TOPO TA cloning 
kit (Invitrogen, USA). The DNA constructs were transferred into 
Escherichia coli TOP10 cells (Invitrogen). The sequencing reaction 
of bacterial 16S rRNA gene clones (20-50 clones per specimen 
and species) and amplified eukaryotic COI, NADH4, 16S, 28S 
and histone 3 genes was performed using a BigDye Terminator v. 
3.1 Cycle Sequencing Kit (Applied Biosystems, USA). Sequencing 
was performed using an ABI PRISM v. 3730 Genetic Analyzer 
(Applied Biosystems). Sequences were proofread using CodonCode 
Aligner v. 3.7.1.1 (CodonCode Corporation, www.codoncode. 
com). Locations of symbionts in gill tissues of Ba. aduloides, Idas 
nsp SAL4 and Idas nsp ESU R were determined using transmis- 
sion electron microscopy (see electronic supplementary material, 
figure SI), following protocols described elsewhere [11]. Symbionts 
were identified by comparing their 16S rRNA sequences with 
GenBank databases using BLAST searches. 

(c) Model selection and reconstruction of the host tree 

DNA sequences from newly obtained hosts were complemented 
with data from GenBank and aligned using Probalign v. 1.4 [12]. 
The best-fitting model of nucleotide substitution was selected for 
each gene and each partition within each gene using the cor- 
rected Akaike information criterion (AICc) in jModelTest 
v. 0.1.1 [13]. We inferred gene trees using the Bayesian 
method implemented in BEAST v. 1.7.2 [14] for each partitioning 
scheme. A Yule speciation model was used as a tree prior. We 
modelled rate variation among lineages using an uncorrelated 
lognormal relaxed clock, with the mean substitution rate fixed 
to 1 to get branch lengths in units of substitution per site [15]. 
Posterior distributions were estimated using Markov chain 
Monte Carlo (MCMC) sampling. Samples were drawn every 
1000 steps over a total of 10 million MCMC steps. Each analysis 
was run four times, with mixing and convergence assessed using 
Tracer v. 1.5. After discarding 10% of the samples as burn-in, 
samples from the four runs were thinned (sampling every 4000 
steps) and pooled together. The best partitioning scheme was 
selected for each gene by comparing marginal likelihoods 
using approximate Bayes factors [16] in Tracer. 

We then performed a combined analysis of all five genes, using 
the selected partitioning scheme and substitution models. Each 
gene was assigned an independent relaxed clock. The protocol of 
the combined analysis was the same as that used for single-gene 
analyses, except the number of MCMC steps was increased to 
20 million. The maximum-clade-credibility tree was drawn from 
the pooled samples. Maximum-likelihood (ML) trees were inferred 
using Treefinder [17], with the same partitioning scheme and 
substitution models as those used in the Bayesian analyses. Boot- 
strapping analyses (1000 replicates) were used to evaluate 
support for the ML tree. 

(d) Molecular dating 

Three fossils were used as calibrations in the Bayesian relaxed- 
clock analysis and were implemented as prior distributions for 
ages of nodes in the tree [18]. 

Mussels from a Middle Eocene (45 Myr) seep deposit in 
Washington State (USA) were recently assigned with some hesi- 
tation to Vulcanidas (as Vulcanidas? goederti Kiel & Amano, 2013). 
This hesitation stemmed from the missing data on anatomy and 
muscle attachment scars, while in all available shell characters 
Vulcanidas? goederti was more similar to Vulcanidas than to any 
other genus of deep-sea mussels [19]. Therefore, we modelled 
the divergence time of the clade including Vulcanidas insolatus 
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Figure 1. Maximum-clade-credibility Bayesian chronograms and estimates of ancestral character states. Pie charts indicate probabilities of each state at nodes 
discussed in the text (yellow rectangles), (a) Chronogram showing the evolution of habitat use. Grey bars are 95% HPD intervals of divergence time estimates. 
Black squares, circles and asterisks at nodes indicate posterior probabilities greater than or equal to 0.99, bootstrap values greater than or equal to 75% (95% for 
lineages highlighted in yellow), and nodes inferred in analyses of both nuclear and mitochondrial genes, respectively. The nodes that were assigned fossils to 
estimate divergence times of the 'childressi' group, and Gigantidas and Vulcanidas genera, are labelled in black circles by numbers 1, 2 and 3, respectively. Abbrevi- 
ations on time scale: Palaeo, Palaeocene; Oligo, Oligocene; P, Pliocene to present; Ma, million years ago. (b) Evolution of the presence/absence of sulfur-oxidizing 
symbionts. (c) Evolution of the presence/absence of methanotrophic symbionts. (d) Evolution of the location of symbionts in the gill epithelium. Squares at tips of 
chronograms b, c and d indicate available data, (e) Evolution of shell lengths, (f) Evolution of depth use. For legibility, and because depth was discretized into nine 
states, a histogram was used instead of a pie chart to display probabilities at the root. Genus names are abbreviated as follows: A, Adipicola; Ba., Bathymodiolius; 
Be., Benthomodiolus; G., Gigantidas', I., Idas; M., Modiolus; 1, Tamu. 



Cosel & Marshall, 2010 from other species (figure la, lineage L9) 
using an exponential prior with a hard minimum bound of 
45 Myr and a mean of 1.2 Myr. This resulted in a soft maximum 
matching the beginning of the Middle Eocene (48.6 Myr). 

Two fossil species from seep deposits on the eastern North Island 
of New Zealand were recently assigned to extant clades based on 
muscle attachment scars. The Middle Miocene (Lillburnian; 12.7- 
15.1 Myr) Gigantidas coseli Saether et al. 2010 was used to calibrate 
the divergence of Gigantidas (figure la, lineage L5) from other 
lineages. An exponential prior with a minimum bound of 15.1 Myr 
and a mean of 1.3 Myr was used for this node, resulting in a soft 



maximum matching the beginning of the Clifdenian stage 
(15.9 Myr). The Early Miocene (Waitakian; 21.7-25.2 Myr) Bathymo- 
diolus heretaunga Saether et al 2010 was assigned to the 'childressi' 
clade (figure la, lineage L4 + L5) and the divergence of this clade 
was assigned an exponential prior with a minimum bound of 
25.2 Myr and a mean of 0.7 Myr. This yielded a soft maximum 
matching the beginning of the Late Oligocene Duntroonian stage 
(27.3 Myr). 

Each of the three fossil calibrations was used in a separate 
analysis, and then we combined them in a single analysis. For 
each analysis, four runs of 20 million MCMC steps, 10% burn-in 



and a 4000-step subsampling of individual runs were used. The 
COI substitution rate and the divergence between Bathymodiolius 
thermophilus Kenk & Wilson, 1985 and Bathymodiolius antarcticus 
Johnson, Won, Harvey & Vrijenhoek 2013, which is thought to 
have occurred at 5 Ma [20], were used to further evaluate the con- 
sistency of results obtained from our three fossil calibrations. 
Molecular dating was also performed using non-parametric rate 
smoothing in r8s v. 1.8 [21,22]. Uncertainties in the fossil calibra- 
tions were modelled using age constraints that corresponded to 
the intervals used in our Bayesian analyses. This analysis was 
run 10 times from different starting points. 



(e) Character evolution and diversification analyses 

Data on habitat use, shell length, depth, presence /absence of 
sulfur-oxidizing and methanotrophic symbionts, and symbiont 
location on the gill epithelium were compiled from published 
sources and added to our new data (see electronic supplementary 
material, table SI). Because depth and depth range vary greatly 
among deep-sea mussels, depth was discretized using intervals 
of 500 m and multiple states were used to account for large 
depth ranges. Transitions between character states were estimated 
using transition matrices for categorical data [23], and a Brownian 
diffusion model was used for natural logarithms of shell lengths 
[24]. These components were added to the analysis that included all 
five genes and fossil calibrations, without the outgroup. As in the 
Bayesian dating analysis, we used four runs of 20 million MCMC 
steps, 10% burn-in and a 4000-step subsampling of individual runs. 

To assess heterogeneity of diversification rates through time, a 
likelihood analysis of speciation and extinction rates [25] was con- 
ducted using the R package LASER. This analysis was based on the 
chronogram estimated by non-parametric rate smoothing of 
the ML tree (free of any speciation model). Distribution of AICc 
differences among a set of speciation models was obtained by 
simulating 5000 trees under the null hypothesis of constant diver- 
sification. The robustness of the results to sampling bias was 
evaluated by running simulations under the assumptions that 0, 
25 and 35% of species were missing. 

A similar simulation approach was carried out using the R 
package TESS to test whether including extinction events in the 
birth- death speciation model improved the likelihood score [26]. 
Simulations were conditioned on the age of the root and the 
number of species in input trees. To assess the impact of methano- 
trophic symbionts and intracellular symbionts on diversification 
rates, character-dependent Bisse [27] and Yule speciation models 
were fitted to the ML tree using the R package diver sitree [28]. 
The significance of likelihood differences between models was 
tested using likelihood ratio tests. These analyses were conducted 
from the ML tree after removing missing data and from the pat- 
terns inferred from Bayesian analyses, assuming sampling biases 
of 0, 25 and 35%. 

Unsupported nodes were collapsed and missing data were 
removed before further analyses. Correlations among categorical 
variables were tested (10 iterations, 1000 simulations) using 
Pagel's method [29] as implemented in Mesquite v. 2.75. Influences 
of habitat use, presence /absence of methanotrophic symbionts, 
and presence/ absence of intracellular symbionts on natural logar- 
ithms of shell lengths and depth were tested using phylogenetic 
analyses of variance [30] with the R package geiger (5000 simu- 
lations) [31]. Correlation between natural logarithm of shell 
lengths and depth was tested using the phylogenetic generalized 
least-squares method [32] with the R package caper. 



3. Results and discussion 

The species tree inferred from the concatenated dataset 
(figure la) revealed support for 10 clades, which were 



consistent with those observed in previous studies [9,33,34]. 
Our Bayesian relaxed-clock analysis, calibrated using fossils 
at three nodes, yielded an estimated mean age of all 
deep-sea symbiotic mussels of 85 Myr (95% HPD interval: 
69-102). Although this result is consistent with their diver- 
gence from other Mytilidae in the Late Mesozoic or Early 
Cenozoic [5,33], there is a lack of mussels at the few methane 
seep deposits and organic falls of this age [35-39]. This dis- 
crepancy may be due to higher levels of homoplasy along 
deep branches and /or to the lack of age calibrations at 
deep nodes in the tree. Alternatively, earliest species might 
be absent from the fossil record [5]. 

Most of the 10 extant clades started diverging from each 
other about 45 Ma (figure la), and the main vent and seep 
clades (L2, L4 and L5) diversified within the last 30 Myr, 
consistent with previous analyses using biogeographic cali- 
bration points [40]. The robustness of these date estimates was 
supported by the consistency of results obtained using each of 
the fossil calibrations separately (see electronic supplementary 
material, figure S2). Differences in mean estimates among 
fossil calibrations were at most 0.38% for the COI substitution 
rate, 0.61 Myr for the divergence between the sister species 
Ba. thermophilus and Ba. antarcticus, 4.05 Myr for the divergence 
of the genus Gigantidas Cosel & Marshall, 2003 (lineage L5), 
and 5.67 Myr for the divergence of the 'childressi' clade (lineages 
L4+L5). In the analysis including the three fossil calibrations, 
the COI substitution rate estimated across the entire tree 
was 1.62 x 10~ 2 substitutions per site per Myr (95% HPD inter- 
val: 1.22 x 10~ 2 -2.09 x 10~ 2 ), a value higher than that estimated 
in vent-endemic annelid taxa, but close to those found in other 
invertebrates from non-chemosynthetic shallow water environ- 
ments [41,42]. Our estimate for the divergence between Ba. 
thermophilus and Ba. antarcticus was 2.77 Myr (95% HPD inter- 
val: 1.7-4.0), a result similar to those obtained from 
population genetic studies [20,43]. 

The inference that the divergences among and within clades 
occurred within the last 45 Myr is remarkable because the eco- 
logical niche of deep-sea mussels — being epifaunal and taking 
up sulfide from the water column — was unoccupied since 
their inferred origin in the Late Cretaceous [44-45]. Even assum- 
ing the lower (younger) bound of the 95% HPD interval as the 
time of origin of the deep-sea symbiotic mussels still implies 
that it took 25 Myr from the origin to the diversification event 
seen in our tree (figure la). This long basal branch led other 
authors to propose the hypothesis that the genus Benthomodiolus 
(L10), which is a sister group to the remaining nine clades, may 
be a relic of lineages that became extinct during the global 
anoxia/ dysoxia event associated with the Palaeocene/ Eocene 
thermal maximum (PETM) around 57 Ma [33,46]. We found 
some support for that hypothesis, since the likelihood of our 
ML and Bayesian maximum-clade-credibility trees under a 
birth- death speciation model significantly improved (p = 
0.002 and p = 0.024, respectively) when including an extinction 
event (5% survival rate) at 57 Ma. However, our simulations 
also show that a scenario without extinction during this time pro- 
duced very similar trees, inducing a lack of statistical power in the 
case of Bayesian trees (figure 2). Overall, our analyses also sugges- 
ted that chemosymbiotic mussels were not very diverse during 
the PETM and that it induced, if anything, a simple slowdown 
of the early diversification rather than a dramatic extinction event. 

The slow initial evolution of deep-sea mussels contrasts 
with the timing of branching events that occurred from 
45 Ma. A likelihood analysis of speciation and extinction 
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Figure 2. Comparison between extinction-free birth -death models (HO) and models including a mass extinction at 57 Ma (H1). Analyses were performed on 
chronograms inferred using (a) non-parametric rate smoothing and (b) Bayesian phylogenetic analysis. Lineages-through-time plots (left) obtained from these 
trees are represented by solid black lines, with additional dashed black lines for the 95% HPD interval estimated from the entire distribution of sampled Bayesian 
trees. Red- and blue-shaded areas correspond to 95% CIs obtained by simulation under HO and H1, respectively, and purple areas show the overlap in confidence 
intervals between models. Distributions of likelihood differences between both models fitted to datasets simulated under HO (red) and H1 (blue) are given on the 
right. Dashed red lines represent HO's 5% rejection levels. Likelihood differences calculated from real trees are given by solid vertical black lines, with an additional 
95% HPD interval for Bayesian trees (vertical black dashed lines). Ma, million years ago; Alog(L), difference between natural logarithms of likelihood. 



rates allowed us to identify a moderate but distinct increase in 
speciation rates, with about 0.17 speciation events per million 
years from 41 to 32 Ma (Middle Eocene to earliest Oligocene; 
figure 3). From 32 Ma to the present day, the estimated 
speciation rate was three times lower. The corresponding 
Yule model with three distinct speciation rates provided a 
significantly better fit (p = 0.007) to the tree than the best 
constant-rate model (i.e. pure birth), even in simulations 
assuming 25% (p = 0.011) and 35% (p = 0.009) sampling bias 
(table 1). Moreover, parts of the gene trees and species tree 
that fall within that interval are consistently characterized 
by poorly resolved nodes and short internode distances 
(figure la). It is worth noting that a complementary analysis 
of speciation rates from 41 Myr to the present suggests that a 
slow decrease of speciation rates, following a Weibull distri- 
bution (/3 = 1.42), might also be an alternative to the simple 
Yule model (x 2 = 8.33, p = 0.004). This indicates that mussels 
have continued to diversify after their initial radiation, as cor- 
roborated by the existence of several young species and 
species complexes [20,46]. 

To test if and how symbiont type and location have 
influenced the diversification of mussels, we estimated the 
ancestral states of these characters. Our results indicate that 
sulfur-oxidizing symbionts are an early acquisition going 
back to the stem of the group, almost 30 Myr before the main 
burst of diversification (figure lb). In contrast, the ability to 
host intracellular symbionts was detected in our tree from the 
Late Oligocene onward, and methanotrophic symbionts only 
from the Middle Miocene onward (figures lc,d and 3). This 
suggests that the evolution of the various symbiotic relation- 
ships did not trigger the Middle Eocene -Early Oligocene 
burst of diversification, although the association with sulfur- 



oxidizing symbionts was probably a prerequisite. A potential 
environmental trigger for this burst might have been the 
appearance of whales, because their carcasses, which produce 
large amounts of sulfide, are thought to be dispersal stepping 
stones [47]. This hypothesis was originally based on the con- 
gruence of molecular age estimates for the chemosymbiotic 
bivalve family Vesicomyidae and the rise of whales, but it 
has been challenged because Eocene and Oligocene whale 
falls lacked associated vesicomyid fossils [48]. However, 
these early whale falls were abundantly colonized by mussels, 
and thus the 'whale stepping stone' hypothesis might apply to 
mussels. Alternatively, the geological record shows a rapid 
increase in the abundance of seep carbonates in the Late 
Eocene [49], indicating a worldwide increase in methane see- 
page. Analogous to the rationale of the 'whale stepping 
stone' hypothesis, this increase in habitat availability and in 
potential dispersal stepping stones could have triggered the 
Late Eocene radiation. A third possible trigger is a drop in 
deep-water temperature beginning in the Middle to Late 
Eocene, associated with the initial glaciation of Antarctica 
[50]. Low ambient water temperatures could have decrea- 
sed the metabolic rate of mussel larvae, thus increasing their 
longevity and enhancing their dispersal capability. This 
would have been an advantage for mussels living in patchy 
deep-sea habitats such as vents, seeps and organic substrates. 

Even if the acquisition of intracellular and methano- 
trophic symbionts was not involved in the Middle Eocene 
to Early Oligocene radiation, symbionts played a role in the 
subsequent evolution of the group. Pagel's tests showed 
that both the presence of methanotrophic symbionts and 
the intracellular location of symbionts are correlated with 
the occurrence at vents and seeps (table 2). This supports the 
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Figure 3. (a) Global and character-dependent lineages-through-time 
plots (log-scaled) estimated from the ML tree smoothed using non-para- 
metric rate smoothing (dashed lines) and from the Bayesian chronogram 
(solid lines). Shaded areas represent 95% HPD intervals estimated from 
the entire distribution of sampled Bayesian trees. Vent and seep lineages 
from figure \a were pooled together in the habitat-dependent plot. 
Colour-coded groups in plots based on the location and the presence/absence 
of methanotrophic symbionts match the distributions of character states 
among taxa in figure \d and 1c, respectively. Vertical dashed lines represent 
shifts of the net diversification rate estimated in the likelihood analysis of 
speciation and extinction rates. 
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Table 2. Tests of correlations between characters. Model used, number of species for which data were available, criteria estimated, degrees of freedom (d.f.) 
and p-values are given for each pair of characters. 



model 


variables 


species 


criterion 


d.f. 


p-value 


Pagel's test 


methanotroph ~ symbiont location 


24 


Alog(L) = 4.916 




0.001 




methanotroph ~ habitat 


36 


Alog(L) = 3.900 




0.024 




symbiont location ~ habitat 


24 


Alog(L) = 4.331 




0.006 


pGLS 


log(size) ~ depth 


48 


intercept: t = 4.460 




0 








slope: t = 0.471 




0.640 








F = 0.222 


46 


0.802 


p An ova 


log(size) ~ habitat 


48 


F = 61.328 


1 


0 








phylogenetic p-value 




0 




log(size) ~ symbiont location 


24 


F = 9.571 


1 


0.005 








phylogenetic p-value 




0.053 




log(size) ~ methanotroph 


36 


f = 3341 


1 


0.057 








phylogenetic p-value 




0.143 




log(depth) ~ symbiont location 


24 


F= 6.994 


1 


0.015 








phylogenetic p-value 




0.097 




log(depth) ~ methanotroph 


36 


F = 6.153 


1 


0.018 








phylogenetic p-value 




0.067 
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adaptive significance of these two character states and the 
hypothesis that symbionts 'increase the metabolic capabilities, 
and therefore the number of ecological niches' colonized by 
mussels [51, p. 475]. We further explored this hypothesis by 
comparing simple Yule speciation models with models in 
which the diversification rates depend on character states. 
Even when assuming high levels of sampling bias, the results 
suggest that diversification rates are higher in lineages host- 
ing intracellular symbionts and methanotrophic symbionts 
(table 3). The acquisition of these two characters is correlated 
with a second wave of mussel diversification at vents and 
seeps (figure 3), and is probably the key to the success of vent 
and seep clades (lineages L2, L4-L6 in figure la). These include 
nearly as many species as from organic falls, but have higher 
speciation rates because they started diversifying only since 
the Miocene. Methanotrophic symbionts were only acquired 
in species hosting intracellular symbionts (figure lc,d', electronic 
supplementary material, table SI) and both characters are 
strongly correlated (table 2). This suggests that the ability to 
host intracellular symbionts was a prerequisite for the 
acquisition of methanotrophic symbionts. 

Our estimation of ancestral habitat types indicates that 
vents and seeps were colonized from organic falls inde- 
pendently in six lineages, consistent with the hypothesis that 
these mussels took 'wooden steps to deep-sea vents' [8]. 
These habitat transitions were unidirectional, suggesting that 
the adaptation to vents and seeps is an irreversible specializ- 
ation. This might be partly explained by the evolution of 
symbioses toward more integrated relationships (i.e. intracellu- 
lar symbioses) and the concomitant reduction of the digestive 
tract in several species [46,52]. Specialization might also 
have resulted from the evolution of other anatomical and 
physiological traits. 

Indeed, we found a significant correlation between habitat 
type and shell size (table 2). We also inferred a trend towards 
larger shells in vent/ seep species from the Early to Middle 



Miocene onwards and towards smaller shells in lineages 
living on organic falls (figure le). This pattern is remarkably 
consistent with the fossil record, which shows that seep- 
inhabiting mussels did not exceed 50 mm in length until the 
Early Miocene [53-55], and then rapidly increased in size 
from 100 mm in the Middle Miocene to more than 300 mm 
today. Compared with the typically small and ephemeral 
organic substrates on the seafloor, vents and seeps provide 
large amounts of energy and habitable space that might explain 
the increased size of species living there [54]. By contrast, 
species living on sparse organic debris would benefit from 
becoming smaller to reduce competition for space and/ or for 
early allocation of resources to reproduction. The latter hypo- 
thesis is corroborated by studies suggesting that mussels from 
organic substrates reach sexual maturity quickly and at much 
smaller size [56-58] than their vent and seep relatives [59]. 

Earlier studies indicated a trend of successive adaptation to 
greater water depth among mussels [33]. This trend is also seen 
in our data for post-Eocene taxa (figure If); however, it 
coincides with habitat transitions from organic substrates to 
vents and seeps, and the latter habitats typically occur in 
deeper water than organic falls (Mann- Whitney: p = 0.001). 
Thus, the colonization of greater depth might just be a conse- 
quence of the habitat transitions. Alternatively, there may be 
a systematic sampling bias because wood and bones are 
much more difficult to locate on the vast deep-sea floor than 
vents and seeps, which can be detected using temperature 
anomalies or seismically reflective sediment layers [60,61]. 

In conclusion, the evolutionary history of chemosymbiotic 
deep-sea mussels shows all the characteristics of an adaptive 
radiation: a burst of diversification into a suite of new habitats 
due to an ecological opportunity, along with physiological and 
morphological adaptations to these habitats. Although these 
mussels rely on their symbionts for nutrition, it is remarkable 
that the initial acquisition of sulfur-oxidizing symbionts 
did not trigger a major radiation. This adds to the growing 
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evidence that an early burst might be rare in adaptive radi- 
ations [62]. The link of the burst of diversification in the 
Middle Eocene to Early Oligocene to ecological opportunities 
such as increased habitat availability or dispersal capability 
raises the question of whether other chemosymbiotic deep- 
sea taxa were similarly affected. Interestingly, a compilation 
and reassessment of molecular age estimates for 14 vent/ 
seep taxa showed that the inferred origins of six of these 
(Alvinellidae, Lepetodrilidae, Alviniconcha / Ifremeria, Provanna, 
Bresiliidae and Bythograeidae) fall within the time interval of 
the mussels' burst of diversification [41]. For the mussels, we 
consider the acquisition of sulfur-oxidizing symbionts as a 
prerequisite for their adaptation to, and successful radiation 
within, chemosynthetic environments. By contrast, the sub- 
sequent acquisition of methanotrophic symbionts allowed 
the colonization of new niches within the vent and seep 
environment, and resulted in a second wave of diversification. 
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